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A simple model of harmonic vibrations in topologically disordered systems, such as glasses and 
supercooled liquids, is studied analytically by extending Euclidean Random Matrix Theory to in- 
clude vector vibrations. Rather generally, it is found that i) the dynamic structure factor shows 
sound- like Brillouin peaks whose longitudinal/transverse character can only be distinguished for 
small transferred momentum, p, ii) the model presents a mechanical instability transition at small 
densities, for which scaling laws are analytically predicted and confirmed numerically, iii) the Bril- 
louin peaks persist deep into the unstable phase, the phase transition being noticeable mostly in 
their linewidth, iv) the Brillouin linewidth scales like p 2 in the stable phase, and like p in the unsta- 
ble one. The analytical results are checked numerically for a simple potential. The main features of 
glassy vibrations previously deduced from scalar ERMT are not substantially altered by these new 
results. 

I. INTRODUCTION 

Recent advances in X-ray and neutron scattering techniques have allowed to obtain very detailed physical insight 
into the high-frequency (0.1-10 THz) vibrational dynamics of supercooled liquids and glasses. Indeed, within this range 
of frequencies their spectra reveal several universal properties [l| , related with the presence of sound- like excitations 
for momenta p of the same order of magnitude of po, the first maximum of the static structure factor (typically 
corresponding to wavenumbers of a few nnr 1 ). In X-ray inelastic scattering experiments, this high-frequency sound 
is mainly revealed as Brillouin-like peaks in the Thz region of the dynamic structure factor, whose position grows 
linearly for p < po, i.e. u)(p) ~ cp (the speed of sound c being quite close to that obtained by acoustic measurements) 
and saturates at a frequency u>(p) ~ too for p ~ Pa- Moreover, the p-dependence of the peak width is often described by 
r(p) = Ap a , where A is basically temperature independent for momenta ranging from O.Olpo to po [3- Interestingly 
enough, T(p) also saturates as the momentum becomes ~ pq. The actual value of the exponent a is still a matter of 
debate among different experimental groups [|| (some proposing a = 2, some a = 4), yet one should bear in mind 
that for somehow smaller momenta, an unambiguous p 2 scaling has been found in optical measurements of sound 
attenuation in amorphous silica 0. 

The fact that the above described features, as well as the Boson peak (see below), are universal, supports the 
hope that most of the underlying relevant physics can be captured by some simple model. Yet, the very nature 
of these systems, intermediate between solids and liquids, poses a considerable challenge to the description of their 
spectra. One could, for instance, take the liquid point of view: the well-known hydrodynamic approximation Q 
predicts the existence of sound-like excitations for all wave numbers, and a sound attenuation coefficient (very much 
the same as the line-width of a Brillouin peak [J) that grows like p 2 . And in fact, as said above, such a scaling 
has been measured for instance in the range 0.01po~0.1j»o in amorphous silica However, when using the value of 
the viscosity of silica at room temperature (an unknown quantity not smaller than 10 13 poises) in the hydrodynamic 
formulas [2; 13, is overestimated by (at least) eight orders of magnitude. Thus, hydrodynamics predicts a completely 
washed-out Brillouin peak, in plain contradiction with experiments. A drastically different approach is to consider 
these excitations (whose inverse frequency is much smaller than the structural relaxation time) as harmonic vibrations 
around a quenched atomic structure, a point of view supported by recent molecular dynamics simulations |2ll I22I l23j ] . 
Given the presence of well formed local structures (Si02 tetrahedra, for instance) a most natural approximation in 
this context is to consider that the oscillation centers form a cr ysta lline structure, the disorder in the atomic positions 
being mimicked by randomness in their interaction potential [jlllO!] (disordered lattice models Q). This approximation 
is particularly appealing for analyzing scattering experiments, since inelastic-scattering from crystals is nowadays a 
well developed discipline. The presence of more than one atomic species in most glass formers produces a complicated 
vibrational structure, with acoustical and optical branches (both longitudinal and transversal) which are degenerate 
in energy for wave numbers close to po- It is clear that a crystalline analogue can be very useful to clean-up the 
mess ^(|. However, disordered lattice models dramatically underestimate the scattering of sound waves 0. 

A somehow intermediate position is held by those studying vibrations around a topolo gica lly disordered Q (liquid 
like) structure. There are basically two such approaches: modified mode-coupling theory | ll| (which is not limited to 
harmonic excitations), and euclidean random matrix theory (ERMT) [T2I IT3L llJ |. ERMT owes its name to the fact 
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that it formulates the vibrational problem as random matrix problem (l5j . The matrices involved are called Euclidean 
random matrices ^(|, and their study has required the development of new analytical tools. Both MCT and ERMT 
predict an enhanced scattering of sound waves as compared to disordered crystals, but up to now had been limited 
(due to technical difficulties) to very simple approximations where the three-dimensional nature of particle vibrations 
was neglected. This paper deals with the extension of ERMT to the physical case of three-dimensional vibrations. 
In this way, we shall learn that indeed the universal features of the high-frequency spectra can be ascribed to the 
topologically disordered structure of supercooled liquids and glasses. However, we must stress that no results are 
available yet for the case were several atomic species are present with very strong local correlations among them (e.g. 
silica tetrahedrons), which would be directly relevant for the interpretation of scattering data. 

In the previous scalar studies, the line width was found to be T(p) ~ Cp 2 + T>p 2 uj 2 , and the peak position uj = cp. So 
for p <C p one gets a = 2, thus recovering the hydrodynamic scaling, but with Brillouin peaks that are still observable 
when p is close to po |l4j . It must be noticed that the possibility remains of a crossover at a momentum p c < po 
(before the broadening saturates) to the a = 4 regime. This possibility should be always considered when discussing 
experimental results obtained under different thermodynamic conditions and for different momentum ranges Q . 

Another quantity accessible to experiments (especially Raman scattering) is the vibrational density of states 
(VDOS), g(w). Its most striking feature is the presence of an excess of states over the Debye lu 2 law in the "low" 
frequency region, i.e. where the dispersion relation is linear (but still in the Thz region). This excess of states is 
seen as a peak when the VDOS is suitably plotted, and has been named Boson peak (BP). There are at least three 
different ways of defining the boson peak from experiments. It is sometimes defined as a peak in Raman scattering 
data, sometimes as a peak in the difference between the observed VDOS of the glass and that of the corresponding- 
crystal (which is cx lu 2 ). The third definition, and the one we will adopt, is to look for a peak in g(uj)/uj 2 . When so 
defined, the peak position to bp usually shifts to lower frequency on heating 01 > except for the case of silica 0|. In 
this material the shift is seen on lowering the density ■ 

Though there are theoretical approaches that explain the BP through anharmonic effects [2(j, there is a growing 
consensus that anharmonicity, although certainly present in real materials, need not be invoked to explain the BP 
or the other features of high frequency sound described above. Indeed, several numerical simulations have shown 
that a model of harmonic vibrations is wholly adequate to describe this frequency range On the other 

hand, even within the harmonic framework the nature of the extra low frequency modes giving rise to the BP is still 
an open point. At a qualitative level, the frequency ujbp is close to the Ioffc-Rcgel frequency loir, the frequency 
at which vibrational modes change from propagating to non-propagating. This is accompanied by a crossover from 
weak to strong scattering of the phonons by the disorder [l(J, |2J| , suggesting the possibility that the excess BP modes 
are localized |25l |26[. However, numerical simulations of amorphous silica have shown that the localization ed ge i s 
at frequencies greater than lo bp and ujjr (27j. This has also been found analytically in the CPA framework |10(. 
What the Ioffe-Regel criterion signals is rather a crossover from phonon-like excitations to a different region where 
the scattering due to disorder is very strong and the modes do not propagate. We call these modes glassons (since 
they do not propagate but "diffuse", they have also been called diffusons |27]|). A large bump of glassons is generally 
found around the Ioffe-Regel frequency, due to the flattenin g of the dispersion relation. This can be considered as the 
glass counterpart of the van Hove singularity of crystals [I(J, l23l l2a | . All the recently proposed theoretical frameworks 
predict that this peak of glassons should move to lower frequencies when approaching an instability transition, where 
negative eigenvalues (imaginary frequencies) appear. This has been related to the arising of the BP [T 
Some recent simulations for silica have stressed that the BP has an strong component of transverse modes |23 

A recent breakthrough of the ERM theory has been to push the analysis of the vibrational spectra of glasses from 
the qualitative level, where mainly orders of magnitude are compared, to a more quantitative one, making sharp 
predictions about the values of universal critical exponents describing the approach to the singularity [30l | . The first 
attempt to confirm the theory by measuring them in a numeric simulation has been quite encouraging |3Cj . The 
emerging picture is the following: the BP modes are given by the hybridization between the phonons and the low- 
energy tail of the glasson peak which softens when the system approaches the instability transition. Such a mechanism 
is strongly suspected to take place in supercooled liquids at the Mode Coupling transition, and there is numerical 
support |31j for the idea that the Mode Coupling transition (actually, cross-over) in supercooled liquids corresponds 
to a smooth phase transition in the energy landscape, from a saddle-dominated region to a minima-dominated one. 

Interestingly enough, since the ERM computation was originally performed in the monoatomic scalar case (i.e. 
the vibrations are all collinear) the results obtained do not depend on the existence of transversal or optic modes 
but describe a general phenomenon occurring when phonons and glassons interact. On the other hand, since in real 
systems the lower energy glassons have been claimed j^, to have transverse polarization, in this work we shall 
support the universal character of the transition by extending the Euclidean Random Matrix Theory to a generic 
model with longitudinal and transverse modes. Our aim is to check the validity of the predictions of the vectorial ERM 
computation on a simple Gaussian model whose spectral properties have been numerically studied by the method 
of moments j^- As a matter of fact the theory predicts that the behaviour of some of the main spectral features, 
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namely the arising of the BP and the broadening of the Brillouin peak, are universal and hence can be captured 
even by the simplest model. We shall show concerning the features mentioned above that the numeric results of the 
Gaussian model agree with the theoretical predictions and are identical to those found in more realistic models of 
fragile glass formers 0, |3(3| . 

The layout of the rest of this paper is as follows. In section II we introduce the theoretical formalism of ERM theory. 
In the next section we discuss the phase transition as it results from an analytical model-independent computation. 
In section IV we check these conclusions by reporting numerical and analytical results of a particular simplified model 
(the Gaussian model). We will see that our analytical results fail in giving the whole shape of the spectrum, due to 
the superposition approximation and to the high density expansion. However, these results describe rather well the 
behaviour of the system near the transition point, with the correct exponents and scaling laws. We summarize our 
conclusions in the last section. 



II. THE EUCLIDEAN RANDOM MATRIX THEORY IN THE VECTOR CASE 



We study a model where particles oscillate around fixed random positions, so that the position of particle i at time 
t is Xi (<) = x^ 9 + ipi (t) ; the x^ 9 are quenched equilibrium positions (whose distribution must be specified) and tfi (t) 
are the displacements. From now on, greek indices will label the Cartesian components of the displacements, <fi(t). 
In the harmonic approximation the Hamiltonian is 



1,N 1,N 1,3 

H [x] = £ v(xi - x,) ~ - £ £ Mmj^tftf 



where the dynamical matrix M is an Euclidean Random Matrix: 



JV 



M^x^-Z^xf -xf) + ^/^(xf-x^), 



(1) 



(2) 



k=l 



with /^(x) = d^ v v(~x). Translational invariancc implies that there arc 3 null eigenvalues corresponding to the rigid 
translations of the system as a whole. 

In the one-excitation approximation (and in the classical limit) the dynamic structure factor measured in inelastic- 
scattering experiments is 



* (1) (p>-) = SE 



5(uj - u>„), 



(3) 



where e n are the eigenvectors of the dynamical matrix and uj n its cigcnfrcqucncies (= square root of eigenvalues). 
The overline means average over the disordered quenched positions, whose distribution P[x eg ] has to be specified. It 
is thus assumed (as it is often the case for disordered systems) that macroscopic observables are self-averaging. The 
density of states (VDOS) is obtained in the limit of large momenta: 



g(uj) = lim 



p— >oo ksTp 2 

We can obtain ^'(p,^ through the resolvent G(p, z): 



N 



e ip-(x° 9 -x^) 


1 




z-M 



Gl(p, z)^r- + Gt(p, z) I 8^ 



PliPv 



(4) 



(5) 



which is an axial tensor that can be separated in a longitudinal term and a transversal one, Gl(p) and Gt(p), 
depending only on the magnitude of p. The dynamic structure factor is obtained from the longitudinal resolvent 
using the distribution identity (x + i0 + ) _1 = P(l/x) — iir5{x) (Plcmclj formula): 



2k B Tp< 



ImG L (pV + iO+) 



(G) 



A transverse dynamic structure factor (not measurable in experiments) can be defined in an analogous way, and will 
have a Brillouin peak corresponding to the transverse excitations. 
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A most important and general result is that for p — > oo the resolvent becomes isotropic: 

1 



M 



= 5^Tt[z-M]-K (7) 



So both longitudinal and transverse structure factors tend to a common limit (the VDOS, see eq. 0} at infinite 
momentum. This implies in principle that both the dispersion relations saturate at the same value. However in 
disordered systems the dispersion relation is ill-defined when u> ~ loir since the Brillouin peak becomes very broad. 
We shall come back to this subject later. 

In order to fully define the model one has to specify the potential w(r), which for the moment we leave unspecified, 
and the probability distribution P[x. eq ]. The choice of P is a crucial one, the various alternatives having a different 
physical meaning. Two cases are especially interesting in the physics of glasses. If one is interested in the instantaneous 
normal mode (INM) spectrum, one should choose 

P[x e? ] ocexp(-^[x e9 ]), (8) 

with (3 = 1/fcgT, since one is studying the vibrations around typical equilibrium configuration |33j . If instead one 
wants to focus on vibrations around stationary points (saddles or minima) of the Hamiltonian with a given energy E, 
the choice would be 

P[^» , E] = j^— S [ Xl - x? (E)} ...S [xjv - x^ (E)} , (9) 

where {x.f(E)} are the Af(E) solutions of the equations <9,iT(x) = under the condition H(ptf) = E. The probability 
function enters in the calculation of integrals of the type 

J ebcidx 2 . . . dyi k g {k) (x : , x 2 , . . . x fe )/(xi, x 2 )/(x 2) x 3 ) . . . /(x fe _i, x fc ), (10) 

where <jrw(xi, x 2 , . . . xj,.) = p~ k fdx.k+1 ■ ■ ■ dxjy P [x] is the fc-particle distribution function. Obviously the computa- 
tions become a lot easier by taking P[x] = 1/V N (i.e. considering an uniform distribution of the equilibriom positions), 
instead of the actual forms JHJ or ©■ This is what we shall do in the following. As we shall later discuss, this brutal 
approximation is not completely useless, as might seem at first sight. For the moment, let us show the main steps of 
the computation in this simplified situation. 

The computation of (J5J is the generalization of the scalar vibrations case 0, so we shall only outline it, 
stressing the differences. Proceeding as in the scalar case, one first expands the propagator in l/z: 

°° 1 1 

GV(P,2) =J2^+ TjvE^ " Xfc) (^ r W. (11) 

r=0 jk 

Consider the term of order r in this expansion. If we do not allow repetitions, there are N\/(N — i — 1)! r* A^ r+1 ways 
of choosing the particle index. There is a factor 1/V r+1 from the average over particle positions, another factor 1/N 
from the definition of the resolvent and an extra factor of V due to translational invariance. Hence terms without 
repetitions contribute 0(p r ) (p = N/V is the particle number density). It is easy to convince oneself that every time 
that we allow a particle repetition, we loose a factor of p. Then, the generic term in the 1 j z expansion is a polynomial 
in p, 

±_Y^e l ^?-< q \Mr) mM = p r A ^(p) + p^A^Hp) + ■■■ + PA^(P). (12) 

jk 

Let us remark that for p = we have 

G MV (0,z) = ^, (13) 



which means that uniform displacements are eigenvectors with zero eigenvalue. This is a consequence of translational 
symmetry, and must hold at arbitrary density. Therefore, we must have A^ s \p = 0) = for all r and s. 

The calculation proceeds by finding the polynomial coefficients A^ s \p) up to a given approximation (1/p expan- 
sion) to all orders r. Then the l/z expansion is resummed. The leading order (that with the highest order in p) 
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defines the bare propagator. The terms with smaller powers of p will take the form of a self-energy, broadening the 
spectral line. Organizing the calculation as explained in ref. ITsl one easily finds A^ r \p) in terms of f^ip), the 
Fourier transform of /^(r), that can be separated in longitudinal transverse parts: 



JUp) = hivf-^ + Mp) (V - ^ 

In the following it will be important to remember that /l(0) = /t(0) = /(0). One finds 



4S r) (p) = ([/(o)-/(p) 



m-hip) 



PfiPu 
p2 



/(o) - Mp) 



8n U - 



PiiPu 



(14) 



(15) 



Thus, at the leading order in the high-density expansion, the 1/z series is geometric and the longitudinal and transverse 
part are independent. The bare propagator (exact in the p — > oo limit) is then 



g^(p,z) = 4%, 2 )^ + 4%, z )(V-^, 

4 0) T (p,z) = [z-p/(0) + p/ L , T (p)]- 1 . 



(16) 
(17) 



Notice that it verifies Ea. i|13|) . as it should. Using © the Brillouin peaks in the dynamic structure factor can be 
recovered in the p — > oo limit. Both the longitudinal and the transverse propagators have a simple pole for each p, 
implying that sound-waves are exact eigenvectors of the dynamical matrix in this limit. The bare propagator describes 
an elastic medium of infinite density where the longitudinal and transverse dispersion relation are given by 

u;l T (p)^( P f(0)~pf L , T (p)) 1/2 , (18) 
which become linear at small p. The density of states, that is obtained from the propagator in the p — > oo limit, is 



rather unnatural for p —> oo |12|: a single Dirac delta function at frequency y pf(0). This pathology is a consequence 

of the lack of a Debye cut-off frequency in our elastic medium (there are infinitely more wave numbers outside an 
sphere of arbitrary large radius, p c , than inside it, but all wave numbers larger than p c will have a frequency arbitrarily 



close to y p/(0)). This problem will be cured by the finite p corrections. 

To calculate the finite p corrections, one must allow for particle label repetitions in the matrix products Ijlip. 
Representing the matrix product by a chain of particle labels, at first order in l/p we have to take care of a single 
particle label repetition, 



1 ... 1 



(19) 



where the unrepeated indices are represented by dots. Calculating the generic contribution of this diagrams for all 
r, one obtains the coefficients A^ r (p), then, resumming the 1/z expansion, one finds (repeated greek indices are 
summed) : 



Gp V (p,z) = G$(p,z) + G$(p,z) 

V^(q, p) = p(/U(q) - / m „(p - q))- 



i 



d 3 q 
p J (2^)3 



y r , A (q,p)Gi°J(q,z)T/ CTp (q,p) 



GfJ(p,z) + o(J-\, (20) 

(21) 



This can be easily interpreted if one has in mind that corrections to the p — > oo limit will take the form of a self-energy: 



1 



«-p/(0)+p/(p)-E(p,«) 



fiiy 



where S(p, z) is a matrix with the standard form 



5V(p, z) = Ei(p)^ + E T (p) 5^ - 

pZ y <- pZ 

Due to (|13fl the self-energy should vanish at p — 0. Now, if the self-energy has a series expansion in l/p, 

E(p,zH£( 1 >(p,z) + £( 2 )(p,z) + ... 



(22) 



(23) 



(24) 
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the contribution of order l/p L being £W(p, z), then the resolvent reads (Dyson resummation) : 

G^(p, z) = G$ (p, z) + G% (p, z)E<JJ (p, (p, *) + G<°> (p, z)X$ (p, z)Gi°j (p, z) + 

GW(p, z)Eg(p, z)Gi°;(p, z)£«(p, z)G(° i )(p, z) + . . . (25) 
Therefore, one is tempted to conclude from eq. lj2"U|) that 

dh 



E#(p,*) = i / ^^A(q,p)Gg(q, 2 r)y (W (q,p) 



(26) 



. . . 1 . 


2 


2 


. 1 . . . 


. . . 1 . 




1 


7, 




1 


1 


..I... 



To see that this is actually the case, one can consider the second order corrections, where we have four different 
topologies: 

(27) 
(28) 
(29) 
(30) 

It is clear that to make further progress one needs some kind of Feynman rules. We have been able to find a set of 
diagrammatic rules valid for planar topologies (i.e. repeated indices can be nested but not intercalated) like eqs. I|27|) 
and (|28(l . which we give below. They should not be used in general situations like ca . (|3()j) . The rules are (see Fig.Q: 

1. Join unrepeated particle indices with a full line. 

2. Join repeated particle indices with a dashed "interaction" line. Momentum is conserved at each interaction 
vertex, so we can attach a momentum to each line. 

3. To a full line carrying momentum k associate a bare propagator Gl°J(k, z). 

4. Place a vertex function V(q, p) at each of the vertices connected by a dashed line. Here p is the incoming 
momentum and q is the one running inside the loop. Beware the vertex function does not commute with the 
bare propagator, so the order they have in the diagram must be respected. The resulting sequence must be 
interpreted as a matrix product. 

5. Integrate over the momenta inside the loops. 

6. Add a factor 1/p for each particle repetition (dashed line). 

In this way, one can easily see that eq. I|27|l represents a term of the Dyson resummation of S^ 1 ', while the other 
three schemes are genuine contributions to S*- 2 '. These rules represent a significant improvement over the previously 
published expansions |l3[ . For example, the repetition scheme (|28[1 involved 16 diagrams in the previous version (see 
appendix of ref. Il3f) . while it now corresponds to just one. The other two diagrams for £( 2 ), which will be neglected 
in the resummation used below (see Fig. |5J) are the diagram corresponding to the particle repetition in l|29|l . 



^2 J 7^3 (|^3 V n<*{P ~ q P) G a/j(*> q)^7(q - k > q)<3!$(£, k)V 5(J (p q, p q + k) x 

xGW(z,p-q + k)V/ TI ,(q-k,p), (31) 

and that corresponding to (|30() , 

1 f d 3 q d 3 k (o) 



P 2 



j (2^3 (^3 Vfla ^ ~ q ' ^' ^ ~ k ' ^ ~ ^ X 



xG^(z,k)V^(p-k,p). (32) 



This three diagrams (see Fig|2J), vanish independently for p = 0, and (in the scalar case) are equivalent to the 39 
diagrams of the previuosly published expansion 

Given these rules, it is a well-known combinatorial result (Dyson equation) that the sum of all planar diagrams 
(cactus approximation) takes the form of a self-consistent integral equation for the self-energy: 



1 f d 3 q 

P J (AT) 



). (33) 
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The first terms leading to this equation are illustrated in Fig. |2J It is most important that this resummation is 
compatible with Ea. l|13[) . An asymptotic model-independent analysis of this equation will be presented in the next 
section, while numerical results for a simple model will be given in section IV. 

Let us remark that the self energy renormalizes the dispersion relations and gives a finite width to the Brillouin 
peaks: 

w l,t(p) = ( w 1t) 2 (?) +R cS l,t(p,^l,t(p)), 

Tl,t(p) = ^^l,t(p,^l,t(p))/ujl,t{p)- (34) 

It is interesting to notice that even if the bare propagator did not have a transverse component the first order 
contribution to the self-energy would generate a transverse excitation band: the scattering among longitudinal 
and transversal phonons is much stronger in our case than in disordered lattices. 

Before proceeding to the analysis of the cactus equation, let us note that the correlations between the equilibrium 
positions of the particles can be taken into account quite easily at the level of the superposition approximation in the 
above approach. This approximation amounts to writing the integral (|10|) as 

jdxidx 2 . . .dx fc 3 (2) (x 1 ,x 2 ).g (2) (x 2 ,x 3 ) . . .3 (2) (x fc _ 1 ,x fc )/(x 1 ,x 2 )/(x 2 ,x 3 ) . . . f(x N -i, x fc ), (35) 

where all the information about the correlations is assumed to be contained pair distribution function 5' 2 '(xi,x 2 ) 
(corresponding to the chosen probability P[x e<? ]). The results derived above for the case without correlations are 
translated to the correlated case (at level of the superposition approximation) by replacing the functions /(x) by 
5 ,( - 2 - ) (x)/(x) [34[. In this way the usual power law divergence of the pair potential for |x| — > is balanced by the 
exponential behaviour of the pair distribution function, and this ensures the existence of the Fourier transform of the 
product /(x)g( 2 )(x). 



III. ANALYTICAL RESULTS: THE PHASE TRANSITION 



In this section we aim to show that from the cactus approximation (equation I33|) it is possible to derive a few 
analytic model-independent results about the arising of the Boson Peak and the broadening of the Brillouin peak. 
These results are expressed in form of scaling laws, whose exponents are predicted in this approximation. As the 
scaling laws obtained do not not depend on the details of the interaction, we claim that they are a general feature 
of topologically disordered systems. It cannot be excluded, however, that the values of the exponents depend on the 
approximation chosen. Simulations and experiments will allow to clarify this point. 

Recalling cq.(0J, the VDOS can be obtained from 

g(u) = - — Im G°°(to 2 + i0+), (36) 

setting in our equation z = uj 2 + i0 + , where G°°{z) = lim^oo G(p, z). 

In order to obtain the quantity G°°(z) one has to solve the integral equation (|33|) in the p — > oo limit: 

= z pf(0) pAG°°{z) - p/t03 / 2 (q)G(q, z) (37) 

where G°°(z), A = (2tt)" 3 / d 3 qf 2 (q) and the last term ( given by the integral of the product of two anisotropic matri- 
ces) in the above equation are matrices proportional to the identity since in the infinite momentum limit everything 
is isotropic. 

In order to deal analytically with (|37|l . the crudest approximation is to neglect the last term, in which case it 
becomes quadratic in G°°(z), and one easily finds a semicircular VDOS, with center at ui = pf(0) and radius 2y/p~A~. 
Clearly enough, this approximation misses completely the low frequency part and it is not suitable to describe the 
phonons. It rather describes the non-propagating but extended modes at higher frequencies (glassons). Though the 
glassons are qualitatively taken into account by this approximation, the semicircular shape of this part of the VDOS 
is an artifact of the approximation, and as we shall see later the complete theory yields a shape more similar to that 
found in real systems. 

By reintroducing in an approximated way the missing term of l|37[) . we are able to study the region of the small 
frequencies. In fact, it is possible to show that it can be written as 
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where b is a positive constant and where the sound velocity c is such that 



1 



1 2 

3Cfji 



CL.T = J- P f'(0)/2 



(39) 



It is not difficult to show that A = a/c 3 , where a does not contain c. Substituting this in eq. (|37[) we obtain a 
quadratic equation for G°°(uj 2 + i0 + ) that gives 



G°°(u 2 +iO + ) ~ (lu 2 - pf(Q) + b/c 3 ) +ipf 2 (0)iu/{4irc 3 ) 



(lu 2 - pf(0) + b/c 3 ) + ip/ 2 (0)w/(47rc 3 



Apa/c 3 



2pa/c 3 



(40) 



In order to study the Boson peak we shall focus in the small frequencies regime. For pedagogical purposes, it will 
very illuminating to consider the case where the VDOS changes because of changes in the density. However, the 
mechanism explaining the arising of the BP is much more general and can be applied even to case where the density 
is kept fixed and other thermodynamic parameters vary. 

Let us consider first the high density limit, i.e. the situation where the first term under the square root is much 
larger than the second. The resolvent becomes 



G°°(u/ + i0+) ~ \w z - pf(0) + -r + i 



,p/ 2 (0) W 



47TC 3 



(41) 



and from eq. I|36(l one gets exactly the Debye's law, g(w) = oj 2 / (2ir 2 pc 3 ) , in the 10 — > limit. As a first step in 
generalizing this result, we note that the same conclusion is obtained by considering the limit of high sound velocity 
(i.e. 1/c 3 small with respect to pf(Q)). 

By decreasing the density we easily show that in eq. f4TJ|l . due to the square root term, G°°(0) develops an imaginary 
part when p < p c , where p c is fixed by (b/c 3 — /O c /(0)) 2 = 4p c a/c 3 . We claim that this can be interpreted as a phase 
transition in the space of the eigenvalues of the Hessian matrix driven by the value of the density. The two phases 
separated from such a phase transition are the stable phase (all positive eigenvalues) and the unstable phase (negative 
and positive eigenvalues). The order parameter is <p = — ImG oo (i0 + ) which vanishes as ip ~ |A|^, with (3 = 1/2. 
The value of the exponent of course could depend on the cactus approximation and is highly remindful of mean field 
theories. Most interestingly, the Boson peak arises, in the stable phase, as a signature of the phase transition. In fact, 
by setting A = (p — p c ) and expanding for small A, 



-Im G°°(uj 2 + i0 + ) oc Imy/aA-iu/u)* 



g(u) ~ uj 3 ' 2 
g(ui) cc 0J 2 /VaA 



if uj/lu* > a A 
if lo/lu* < a A 



(42) 



where a is a positive constant and lo* = 2ir \J aj p c c 3 . Hence there exists a frequency that signals a crossover from a 
Debye behaviour to a different kind, namely g(uj) ~ lo 1 , 7 = 3/2. We identify the Boson peak with such frequency, 
i.e. ll>bp = olAlo* . 

From the experimental point of view, this implies that the BP is indicated from a peak in the function g(uj)/uj 2 m 
not in g(uj). The BP frequency moves toward when approaching the transition (from the stable side) and its height 
diverges. Eq. Q shows that at the level of the one-phonon approximation it can also be detected in the p — ► 00 limit 
of the dynamic structure factor S(j>, uS). 

As said above, the control parameter need not be the density. For example, by fixing the density and letting the 
sound speed to vary, we find that 



ip oc Im 



^{b-pc 3 I{0)f 



4pac 3 



(43) 



Then the critical sound speed is fixed by the condition (6 — pc 3 /(0)) 2 = 4pac 3 . Now A = c — c c and once again 
ip = a'\A\P with /3 = 1/2. 

Since the quantities A, b and /(0) depend on the thermodynamic parameters, in principle each of them could be 
chosen as a control parameter to describe this phase transition. However, one parameter (the potential energy) seems 
to be more physically meaningful than the others when the energy landscape approach to the glass transition is taken 
into account. It has been shown indeed in numerical simulations that the typical stationary point of the Hamiltonian 
closest to equilibrium configurations is a saddle above the mode-coupling temperature and a minimum below |3lj| . 
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Hence eq. (|37|) describes (in an approximate way) this phase transition in the space of all stationary points (saddle 
and minima). This we shall call the saddle-phonon transition. 

Since the number of eigenvalues of a stationary point depends only on its energy (in the thermodynamic limit) |3l| . 
it plays a special role as a control parameter in the study of the Boson peak. Stationary points (and their energy), 
though, are easily accessible in simulations but not in experiments. It is thus better to formulate the scaling laws 
arising from the theory without any reference to the parameter chosen. We then recast eq. (|42|l as 

^A)=^M-A-), *(*)~{^. HI , (44) 

with A defined in terms of an arbitrary control parameter. Since the Boson Peak signals a crossover between the two 
phases, this scaling law implies that uj^p ~ A p and g(ujgp, A)/u> 2 3P ~ A _7? , with 

V = P(2 - 7) • (45) 

Since the result is model-independent this law is expected to hold very generally, and in the next section we will check 
that numerically in a simplified model. The ERMT (in the cactus approximation) predicts p = 1, 7 = 3/2, rj = 1/2. 
Numerical results on a fragile glass former have already turned out to be in a reasonable agreement with such scaling 
laws |3C| . 

A similar asymptotic analysis gives the scaling behaviour of the linewidth T(p) of Brillouin peaks. Interestingly 
enough, it changes at the saddle-phonon transition. From the analytic point of view, the longitudinal and transverse 
T(p) are given by the corresponding imaginary part of the self-energy computed at the position loq of the peak. In 
both cases the leading contribution in the limit of small momenta is obtained by considering the large q contribution 
in the integral of eq. (|33|1 : 

Im S (p,^-^/^-^ (q ,p). (46) 



2lo J (2tt) 3 

In the stable phase the Debye regime (g(uj) ~ uj 2 ) at low enough frequencies (uj < lobp) implies 

T L , T (p)~g(Lu)p 2 /Lu 2 ~p 2 . (47) 

This is the asymptotic result both in the longitudinal and the transversal case, but in practice one should be very 
careful when fitting the experimental results with the above law. In fact, when approaching the transition from the 
stable phase, the Boson peak frequency u>bp shifts to zero and the Brillouin peak might fall in the region where 
g(ui) oc uj 3 / 2 yielding a flattening of the scaling law (the exponent would be 3/2 instead of 2). Moreover, a crossover 
from p 2 to a p 4 behaviour in the region where the dispersion relation is still linear cannot be ruled out. In general, 
both the p 2 and p 4 terms are present, and their relative weights depend on the thermodynamic parameters flsTITij. 

In the unstable phase instead one has g(uj) ~ uj, due to the fact that ImG°°(0) ^ and to eq. (|3(j|) . then the 
broadening becomes 

r L , T ~ p. (48) 

This theoretical result is very suggestive because the possibility arises of investigating the saddle-phonon transition 
through measurements of the Brillouin peak. 

We stress that the results of this section are model-independent. Next section we will test them in the simplest 
case, the Gaussian model. 



IV. NUMERICAL RESULTS FOR A GAUSSIAN MODEL 



Here we solve numerically the cactus equation for the case where f(p) has a Gaussian form and compare with 
direct numerical results for the same model. This will confirm that the saddlc-phononn transition described by the 
Euclidean Random Matrix theory is not an artifact of the approximation involved (cactus resummation). The model 
is described by 

iUp) = h(vP^ + hip) (<W - P ~f^j , Il,t{p) = Q) 3/2 cxp (-p 2 /2al iT ). (49) 

This choice for f(p) is mainly due to its simplicity. However, we have shown that many features, ranging from the 
behaviour of the Boson peak close to the saddle-phonon transition to the width of the Brillouin peak in the stable 
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and unstable phases are independent of the details of the model. Hence we expect the Gaussian model to be as good 
as any other to study such features. Moreover, as we discussed at the end of sec. II, the superposition approximation 
amounts to taking favip) = F\9{ r ) v nv( r )\, meaning that /l,t{p) are finite at p = 0. Thus a Gaussian approximation 
for f(p) at low enough momenta is always possible. As an example, in Fig. (top) we show the longitudinal and 
transverse components of the Fourier transform of g{r)Vfi U {r) for a soft sphere system (i.e. a one-component liquid 
with pair potential v(r) = r~ 12 ) and the comparison with a Gaussian fit. The fit gives <jt — 2<7z, (hereafter we shall fix 
<tl to Co, the lenght scale of the problem). We can check that the Gaussian model is able to reproduce the qualitative 
behaviour of a real system by looking at Fig. where it is shown that the density of states of a Gaussian euclidean 
random matrix looks qualitatively like the INM spectrum of a liquid. Thus a Gaussian model is not an outrageous 
approximation. 

We shall consider various values of the density, which is here a possible control parameter, comparing the analytical 
(cactus) results with the numerical spectra and dynamic structure factor obtained from the method of moments |32j | . 
This method allows one to obtain the density of states and the dynamic structure factor of a given N x N matrix up to 
TV = O(10 6 ) in a reasonable CPU time. The numerical solution of eq. I|33|) . which is actually two (coupled) equations 
for the longitudinal and transverse parts of the resolvent, requires some work. It can be solved by iteration at fixed z 
as long as one can do the three-dimensional integral sufficiently fast. To do so, we write the equations as convolutions 
by expanding the vertex V(q, p): this gives four terms, which can be evaluated with two convolutions. These can 
be calculated very efficiently using the fast Fourier transform (FFT), as long as the functions are sufficiently smooth. 
Since this is not the case, the resolvent must be separated in a regular part which tends to zero at infinity, and a 
part which is quasi-singular and has a finite limit. After the convolution of this part is worked out analytically, the 
remaining convolutions can be done using FFT. In this way, each iteration takes a time 0(Alog A), and an iterative 
scheme can be reasonably implemented. 

A. p> p c : the stable phase (phonons) 

In the high density regime we are in the phonon or stable phase, since all the eigenvalues are positive. Fig. [S] 
shows the VDOS, and its Debye behaviour (oc u> 2 ) for small frequencies. In this regime the approximations used in 
deriving the integral equation (|37[) are quite good since the analytic solution reproduces the numerical spectrum (and 
in particular the Debye behaviour) rather accurately. The large peak at high frequency arises from a pileup of states 
due to the flattening of the dispersion relations. Those excitations are the non-propagating glassons described by the 
high-frequency limit of cq. i|37|) . In a way, the peak can be considered as the off-lattice analogous of the van Hove 
singularity, the main difference being that density fluctuations smooth the cusp, giving it rather a semicircular form. 
Finally, let us note that the analytical solution misses the high-frequency tail of the VDOS. These modes are expected 
to appear in a number which is exponential in p, and thus cannot be recovered in a l/p expansion. 

In Fig. we plot the longitudinal and transverse components of the dynamic structure factor for three values of 
p. The agreement between numerical and analytical solutions is better at the higher values of p, though the position 
of the Brillouin peak is correctly reproduced at all p. Since the peaks are well defined, we are allowed to study 
the dispersion relations as well as the linewidth T as a function of p. They are obtained by fitting the peaks with 
,uj) oc uj\ t (jp)Tl^t{p)I [(^ 2 — uj2 l,t{p))' 2 + cj2 r|,T(-P)] > following the experimental procedure Q. While the 
agreement about the numeric and analytic dispersion relations is striking, there remains some discrepancy in the 
linewidth, especially at the lower momenta of the longitudinal case. It is worthwhile to note that both dispersion 
relations saturate at the same value at vary large momenta, though they can be very different in the regime which is 
presently explored by simulations and experiments (pa® ~ 1). Interestingly enough, both relations coincide with the 
bare dispersion relation ujq(p) of ea. (|18fl . the renormalization due to the self-energy in l|34[) seeming negligible. As for 
the linewidth, the numeric and analytic (integral equation) results are in reasonable agreement with the asymptotic 
theoretical prediction r ~ p 2 . Furthermore, the scaling seems to be independent from the longitudinal or transverse 
nature of the excitations. However, it cannot be ruled out from our data the possibility that a crossover between p 2 
and p A is verified before saturating at the limiting value. 

B. p ~ p c : near the transition. 

In Fig. [S] we show the spectrum of Gaussian matrices (in terms of eigenvalues as well as frequencies) as obtained 
from the method of moments for several values of the density near the critical one, which turns out to be p c ~ 0.54. 
If we lower the density we enter the unstable region and this is revealed by the appearance of an extensive fraction of 
unstable modes (imaginary frequencies). A BP appears near the transition point in the low frequency region. 
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We now obtain the exponents of the transition (Fig.[5J). Figs.Ek andEJ; show that the position of the BP is linear 
with respect to A = (p — p c ) and that the height of BP diverges as A -1 / 2 . This confirms the theoretical predictions 
v = 1 and T] = 1/2. In Fig. Eh we determine the value of 7 by studying the fraction of unstable modes. Indeed, from 
eq. in the region of parameters such that p < p c , we can argue that the VDOS behaves as g(u>) = u 1 g(u/\ A|), 

where g(x) is a scaling function and 7 = 3/2. The fraction of unstable modes is defined as /„ = f_ g\(X)d\. We 
thus have that 

poo 

f u (A)= du>u>ig(w/\A\)~ |A| 1+ ^. (50) 
Jo 

We find f u ~ (p c — p) 5 / 2 , i.e. 7 = 3/2. Finally, the order parameter ip vanishes as {p c — p) 13 with (3 = 1/2 (Fig.^Ji). 

The phase transition is also found by the numerical solution of the integral equation (eq. I33J) . as one can see 
from Fig. I1UI Actually, this solution does not give the spectrum very accurately (one can still discern the semicircle 
centered in to ~ p/(0) rather separated from the rest of the spectrum, while we have seen that they mix at low enough 
density — Fig.[BJ and it overestimates the critical density, p^naiyt ~ 1.98 > p" umcnc . In Fig.^Jwe check that the order 
parameter vanishes with the correct exponent, j3 = 1/2, and also for the fraction of unstable modes. It is interesting to 
remark that these exponents are also found for a larger class of nearly Gaussian matrices, corresponding to functions 
f[p) oc exp (-p 2 {l + ap 2n )/2a 2 ). 

C. p < p c : the unstable phase (saddles). 

From the results of the integral equations, we show in Figs. 1121 and 1131 that the Brillouin peaks still exist in the 
unstable phase (below p c ). The main difference with the stable phase is the peak around to ~ 0, due to the negative 
eigenvalues. Of course, since ERMT is a purely harmonic theory, it is expected to describe correctly only frequencies 
much greater than the inverse of the a (structural) relaxation time r Q . So as long as wr Q 3> 1, our computation should 
be reliable to study also the saddle phase. In fact, our finding of Brillouin peaks even when negative eigenvalues are 
present is consistent with the fact that in real systems Brillouin peaks related to the high-frequency sound are found 
even well above the mode-coupling temperature p], which corresponds to the saddle-phonon transition. 

Furthermore it is worthwhile to notice that, as predicted from the analytic asymptotic analysis, the saddle-phonon 
transition produces a change in the linewidth scaling. In Fig. 1131 we show that from the integral equation one gets 
r(p) ~ p. It would be very important to confirm this result in simulations of realistic models or, better, in experiments. 
Finally, at fixed momentum, the scaling laws ujq ~ p 1 ^ 2 for the position of the peak, and T ~ p^ 1 for its width remain 
unchanged across the transition (Fig. I12f) . 

D. The sound velocity as a control parameter 

Finally, we look at the instability transition as driven by the speed of sound, as discussed at the end of section III. 
Since c oc (1 + 2(cl/ct)~ 3 )~ 1 ^ 3 , this is determined by the ratio between longitudinal and transverse sound velocities, 
i.e. by the anisotropy. We thus have studied our ensemble of Gaussian matrices at a fixed density p = 1 as a function 
of the ratio cl/ct- The control parameter is now A = (c — c c ), and we find that c c — 0.58. In Fig. ^]we see that 
a BP appears on approaching the critical value, which signals the transition to a region with an extensive number 
of negative eigenvalues. Fig. 1151 shows that the critical exponents are the same as when the density is the control 
parameter. This confirms the universal character of the phase transition discussed before. 

V. CONCLUSIONS 

In this work we have extended ERM theory to the physical case in which both longitudinal and transverse vibrations 
are present. Even if the model is exceedingly simple (harmonic vibrations around fully disordered oscillations centers) 
it has many features in common with real glasses and supercooled liquids: 

1. A stable phase (no negative eigenvalues) is present, where low energy excitations can be considered as phonons, 
in the sense that they have a well defined wavenumber and that a Debye density of states is found (g(u>) oc lo 2 ). 

2. There is a phase transition from the stable phase to a saddle-phase, which in our model can be controlled by the 
density. Numerical simulations |3l| of supercooled liquids have strongly suggested that such a phase transition 
occurs at the mode-coupling temperature. In our model, the transition leaves a clear mark in the low energy 
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part of the spectrum in the stable (minima) phase: an excess of states, that can be described with scaling laws 
and critical exponents. We propose to identify this spectral feature with the Boson peak of supercooled liquids 
and glasses. This implies that our scaling laws should apply to experimental spectra. We have also found scaling 
laws for the approach from the unstable (liquid) side of the transition, but it is still not clear how they could 
be confronted with experiments. 

3. In this stable phase, the scattering of sound waves is stronger than in disordered crystalline systems: the 
linewidth scales like p 2 even in the harmonic approximation, and there is a significant mixing among longitudinal 
and transversal excitations. It is interesting that, according to our calculation, the scattering of sound waves 
is stronger in the unstable phase. To be precise, the (harmonic) linewidth is found to be T(j>) ~ g(u)p 2 /u> 2 , so 
that in a liquid phase, where g(tu) ~ ui, one has T(p) oc p rather than p 2 . 

We have studied our model using a high-density expansion, and we have compared it with extensive numerical 
simulations. As one could expect, the comparison is rather good at high density, in the phase where the system is 
stable. On the numerical side, we have been able to study in detail the vibrational spectrum in the neighborhood of 
the phonon-saddle transition. Indeed, previous liquid simulations |30t l3lj had been carried out in very small systems 
where low- frequency, long-wavelenght excitations were not present, and our previously studied scalar model did not 
even have a phonon-saddle transition. The present approximation overestimates the critical value of the density by 
a factor four, a tremendous improvement over the calculation for the scalar case, where an inexistent transition was 
found due to the approximations. Since we arc dealing with a high density approximation it is not surprising to find 
significant differences with the numerical calculations at low densities. Nevertheless, universal features like scaling 
laws and critical exponents are captured correctly, as the comparison with the numerical simulation shows (surpringly, 
even the oversimplified scalar model did capture the value of exponents). In any case, considering transverse modes is 
a necessary step towards a quantitative theory of glass vibrations. In fact, in systems with purely repulsive potentials 
(like our Gaussian model, or soft spheres) it is vector displacements that permit the rising of an instability, which 
would be otherwise absent. For the Gaussian model, wc have indeed verified that the instability transition can also 
occur controlled by the ratio of the sound velocities, and that the scaling laws apply as well with this control parameter. 

We think that, from the theoretical side, the subject is mature for a detailed experimental study, since one has 
scaling laws and critical exponents that will hopefully describe the experimental vibrational spectra. Yet there are 
important experimental difficulties related with the fact that the vibrational spectra could be reasonably described as 
harmonic only at low temperatures. This is best explained in terms of the potential energy landscape. A numerical 
simulation of a supercooled liquid typically gets out of equilibrium at the mode-coupling temperature. This means 
that the system is exploring the energy minima close to the phonon-saddle transition where the Boson peak should be 
a very prominent feature of the spectrum. Unfortunately, in experiments one cannot cool as fast as in a simulation, 
meaning that the system goes out of equilibrium rather at the glass temperature (roughly speaking, the system goes 
too far into the stable phase) . Instead, what we propose is to study the harmonic spectrum (experiments done at low 
temperatures), of systems that have fallen out of equilibrium at temperatures near the mode-coupling temperature. 
We expect that this can be achieved by using the new ultrafast cooling techniques (hyperquenching) . Preliminary 
results |36l | indicate that indeed the Boson peak is enhanced at faster cooling, in qualitative agreement with our 
expectations. However, much work is needed to check to what extent the present results can quantitatively describe 
the experimental spectra. In particular, the precise value of the critical exponents could depend on the algebraic or 
exponential decay at large distance of the interaction potential. This is not the case in the present approximation. 

Finally, an important result of the vectorial analysis is the fact that the resolvent becomes isotropic in the infinite 
momentum limit. This means that both dispersion relations saturate at the same value (corresponding to the glasson 
region of the VDOS), and that both structure factors tend to the VDOS for p — > oo. This, strictly speaking, precludes 
the identification of BP modes as those belonging to the end of the transverse branch |2^|. However, since the 
dispersion relations can reach the limiting value in very different ways, for intermediate values of p (those presently 
studied in simulations) they can be very different and may seem to saturate at different values. As an example, in 
Fig. |21 (bottom) we plot the bare longitudinal and transverse dispersion relations of a monoatomic soft-sphere liquid. 
Assume (as is the case for the Gaussian model) that they arc not significantly altered by higher order terms. Then 
we see that the longitudinal branch shows larger fluctuations around the limiting value than the transverse one and 
that they become similar only at very large momenta. But at these valus of p scattering between the two branches is 
so large that the question of the longitudinal or transversal nature of a given mode is rather ill-posed. In any case, to 
better assess the role of the different modes in more complex material, a model including optical branches is needed. 
Thus we plan to investigate a binary model in the near future, as well as to consider more realistic interactions within 
the same formalism, along the lines discussed at the end of sect. II. 
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p/po 



FIG. 3: Top: The Fourier transform (longitudinal and transverse components) of p(r)9 M ^w(r) (g(r) is the radial distribution 
function and 9 M „u(r) is the second derivative of the pair potential) as obtained from a simulation of a soft sphere system at 
T — 16. The Gaussian fits (also shown as solid lines) are made in the small momentum region, and the ratio of two variances 
(which is in our model the square of the ratio of the longitudinal and transverse sound speed) is ~ 2. This is the value we fixed 
in the model. Bottom: The bare dispersion relations (see text) oj° t(p) = p[/(0) — fip)] 1 ^ 2 corresponding to this system. They 
saturate at the same value but the asymptotic behaviour is reached for values of the momenta greater than the ones studied in 
experiments. 
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FIG. 4: Top: The density of states for our ERM model in the unstable phase (p = 0.2), as obtained numerically from the 
method of moments. This numerical method allows to build the density of states of a large matrix (here the Hessian of a 
N = 10 6 particle system) by a sum of delta-functions. The width of these is responsible for the oscillations in the resulting 
VDOS that should vanish once one averages over a big number of samples, ujq is the saturation frequency for high values 
of momenta in the dispersion relation. Bottom: The INM spectrum as arising from a simulation of a soft sphere system of 
N = 2048 particles at T = 0.68, i.e. in the liquid phase. We average here over O(10 2 ) samples, so there are not big oscillations 
in the VDOS. As before, ui is the frequency corresponding to the saturation of the dispersion relation. In both cases the 
density of states for imaginary frequencies is plotted in the negative real axis (ico — > — o>), as usual. 
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FIG. 5: The VDOS g(w) as a function of eigenfrequencies divided by the predicted Debye behaviour u} 2 for p = 4 > p c . In the 
inset, we show g(co) vs. u> both numerical (obtained via the method of moments) and analytical (from the numerical solution 
of the cactus equation 




FIG. 6: The longitudinal (top) and transverse dynamic structure factor in the high density (stable) region for three values of 
exchanged momenta. The comparison between analytical (solid line) and numerical (dashed line) data shows that there is some 
difference especially for small values of momenta. The range of frequencies is in the Debye regime (see Fig. ■ 
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FIG. 7: The dispersion relations and the linewidth of the Brillouin peaks for p — 4 ^> p c , both from numerical evaluation and 
from analytical computation feci. I33P . 




FIG. 8: (a): The spectrum of the Gaussian ERM as obtained numerically with method of moments for several values of the 
density p, below and above the critical density p c ~ 0.54. All values of the densities are given in units of cr^" 3 (see eq. (1491 ). 
(b): Same data, a zoom of the interesting region near A = 0. The value of g(0) grows on lowering the density, (c): Plot of 
g(u>)/uj 2 vs. lo for densities p > p c . The Boson peak is defined as a peak in this plot; we see that a BP appears on lowering the 
density, (d): The VDOS as a function of real (positive axis) and imaginary frequencies (negative axis). 
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FIG. 9: Data taken from numerics (method of moments). The critical density in the following fits has been fixed to p c = 0.54 
and capital letters are the fitting parameters, (a): The position of Boson peak is plotted as a function of the density near the 
critical point. We find that lobp vanishes linearly in A = p — p c . (b): The fraction of unstable modes vanishes as (p c — p) 2 ' 5 , 
thus giving a value 7 = 3/2. (c): The height of the BP, defined by g{u>Bp)/u BP , diverges as A~ v , with 77 = 1/2. (d): The 
order parameter ip = — ImG°°(0) (see text) vanishes as (p c — p) 13 ', with (3 = 1/2. 
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FIG. 10: Top: The density of eigenvalues obtained from the numerical solution of equation l-i'-il at fixed cl /ct = 2 for several 
values of the density p. The transition from a stable phase to an unstable phase is clear since an extensive number of negative 
eigenvalues appears. Bottom: The transition can be seen also by plotting the spectra as a function of the frequencies, with the 
unstable modes (i.e. imaginary frequencies) on the negative real axis. 
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FIG. 11: The order parameter if = -Im G°°(E=0) fhere obtained from the numerical solution of the self-consistent equation (1331) 'l 
behaves as (p — Pc) 13 , with (3 = 1/2, as predicted by the asymptotic analysis of the equation. In the inset, we show that the 
fraction of unstable modes vanishes at the transition as a power law with the exponent 7 + 1 = 5/2. 




FIG. 12: The position of the Brillouin peak loo, here obtained from the numerical solution of the integral equation 1331 scales 
as y/p in a wide range of densities across the critical one. The linewidth F decreases as 1/p in the same range of densities. 
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FIG. 13: The structure factor in the unstable phase, from the solution of eq. (1331 . shows Brillouin peaks (left); the dispersion 
relation (top right) obtained from these peaks is similar to that in the stable phase. The linewidht scaling, however, is changed 
to T ~ p (bottom right). 
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FIG. 14: Data obtained from the solution of eg. 1331 At fixed density p — 1 , the VDOS (evaluated from the integral equation) 
depends on the ratio between longitudinal and transverse sound velocity in the low energy region. This is exactly the same as 
varying the value of the sound speed, since c tx (1 + 2(cl/ct)~ 3 )~ 1 ^ 3 ■ 



22 



0.040 



0.035 



0.030 



9- 0.025 



0.020 



0.015 



0.010 



0.06 



0.05 
0.04 
0.03 
0.02 
0.01 
0.00 



~i 1 1 1 1 r 



••• 



0.50 0.51 0.52 0.53 0.54 0.55 0.56 0.57 
C 



c =0.577 



0.07 



0.08 



0.09 



0.1 



(c-c) 



1/2 



0.11 



0.12 



0.13 



0.14 



FIG. 15: At fixed density (p = 1 here) the transition is driven by the speed of sound, and the order parameter vanishes as a 
square root. In the inset we show that the fraction of unstable modes vanishes as (c — c c ) 2 ' 5 , as predicted by ERMT. The data 
are obtained from the numerical solution of the cactus equation. 



